Let’s look at the ℓ = 0 spherical Bessel function and find their ``zeros.’’

Plot[SphericalBesselJ[0,r],{r,0,15},PlotRange->All]

The ``zero’’ points  for ℓ = 0 are easily found to be, π, 2π, 3π, etc.   for  𝓃 = 1, 2, and 3,  etc.

Let’s look at the ℓ = 1 spherical Bessel function and find their ``zeros.’’

Plot[SphericalBesselJ[1,r],{r,0,15}]

Find the first 3 zeros of (Subscript[𝒿, 1])((Subscript[𝓀, 𝓃])𝓇) for 𝓃 =1, 2, and 3

r1=Table[FindRoot[SphericalBesselJ[1,r]==0,{r,(2 n+1)π/2-0.1}],{n,1,3}]

m=Table[r1[[p,1,2]],{p,1,Length[r1]}]

Draw the spherical Bessel functions for  ℓ=1 , 𝓃 = \[ScriptOne], \[ScriptTwo], and \[ScriptThree] 
and fit them into region  \[ScriptZero] < (Subscript[𝓀, 𝓃])𝓇<\[ScriptOne]

list1=Table[SphericalBesselJ[1,r m[[i]]],{i,1,Length[r1]}]

Plot[list1,{r,0,1}]

Draw the first 3 wave functions for  ℓ=\[ScriptOne] , 𝓃 = \[ScriptOne], \[ScriptTwo], and \[ScriptThree] 
that fit in region  \[ScriptZero] < (Subscript[𝓀, 𝓃])𝓇<\[ScriptOne]        (Subscript[𝓊, 𝓃ℓ])(r) = (Subscript[A, 𝓃ℓ]) 𝓇 (Subscript[R, 𝓃ℓ])(r)  = (Subscript[A, 𝓃\[ScriptOne]]) 𝓇 (Subscript[R, 𝓃\[ScriptOne]])(r)

list2=Table[r SphericalBesselJ[1,r m[[i]]],{i,1,Length[r1]}]

Plot[list2,{r,0,1}]

Let’s look at the ℓ = 2 state spherical Bessel function and find their ``zeros.’’

Plot[SphericalBesselJ[2,r],{r,0,15},PlotRange->All]

r2=Table[FindRoot[SphericalBesselJ[2,r]==0,{r,(n+1)3}],{n,1,3}]

m2=Table[r2[[p,1,2]],{p,1,Length[r1]}]

Draw the spherical Bessel functions for  ℓ=2 , 𝓃 = \[ScriptOne], \[ScriptTwo], and \[ScriptThree] 
that fit them into region  \[ScriptZero] < (Subscript[𝓀, 𝓃])𝓇<\[ScriptOne]

list3=Table[SphericalBesselJ[2,r m2[[i]]],{i,1,Length[r2]}]

Plot[list3,{r,0,1}]

Draw the first 3 wave functions for  ℓ=2 , 𝓃 = \[ScriptOne], \[ScriptTwo], and \[ScriptThree] 
 in region  \[ScriptZero] < (Subscript[𝓀, 𝓃])𝓇<\[ScriptOne]        (Subscript[𝓊, 𝓃ℓ])(r) = (Subscript[A, 𝓃ℓ]) 𝓇 (Subscript[R, 𝓃ℓ])(r)  =  (Subscript[A, 𝓃\[ScriptTwo]]) 𝓇 (Subscript[R, 𝓃\[ScriptTwo]])(r)

list4=Table[r SphericalBesselJ[2,r m2[[i]]],{i,1,Length[r2]}]

Plot[list4,{r,0,1}]

Normalize the 3 wave functions (Subscript[A, 𝓃\[ScriptTwo]]) r (Subscript[𝒿, \[ScriptTwo]])((Subscript[𝓀, 𝓃])r) for 𝓃 = 1, 2, and 3

norm2=Table[Solve[(∫_0^1)((list4[[i]])^2) (norm^2)\[DifferentialD]r==1,norm],{i,1,Length[r2]}]

norm2[[1,2,1,2]]

list5=Table[norm2[[i,2,1,2]]*list4[[i]],{i,1,3}]

Plot[list5,{r,0,1}]